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Modification of the 21-cm power spectrum by quasars 
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ABSTRACT 

We assess the effect of a population of high-redshift quasars on the 21-cin power spec- 
trum during the epoch of reionisation. Our approach is to implement a semi-numerical 
scheme to calculate the three-dimensional structure of ionised regions surrounding 
massive halos at high redshift. We include the ionising influence of luminous quasars 
by populating a simulated overdensity field with quasars using a Monte Carlo Markov 
Chain algorithm. We find that quasars modify both the amplitude and shape of the 
power spectrum at a level which is of the same order as the fractional contribution 
to reionisation. The modification is found both at constant redshift and at constant 
global neutral fraction, and arises because ionising photons produced by quasars are 
biased relative to the density field at a level that is higher than steller ionising photons. 
Our results imply that quasar ionisation will need to be included in detailed modelling 
of observed 21-cm power spectra. 
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1 INTRODUCTION 

The reionisation of cosmic hydrogen is commonly believed 
to have been due to UV photons produced by the first stars 
and quasars and was an important milestone in the history 
of the Universe (e.g., iBarkana fc Loebll2001^ . Reionisation 
starts with isolated regions of ionised hydrogen (H II ) form- 
ing around galaxies and quasars, which later grow and merge 
to surround clusters of galaxies. The reionisation process is 
complete when these H II regions overlap and fill the volume 
between galaxies. Various experiments are currently under- 
way to measure 21-cm emission from the pre-overlap inter- 
galactic medium (IGM) and thus observe the evolution of 
the ionisation stucture directly. These experiments include 
the Low Frequency Arraj0 [LOFAR], the Murchison Wide- 
field Arra^Q (MWA) and the Precision Array to Probe Epoch 
of ReionizatioiJl (PAPER). 

Several probes of the reionisation epoch in redshifted 
21-cm emission have been suggested, which include: obser- 
vation of emission as a function of redshift averaged over 
a large area of sky; imaging of individual Hll regions; ob- 
servation of the power spectrum of fluctuations. Observing 
emission as a function of redshift averaged over a large area 



* Email: pmgeil@unimelb.edu.au 
^ http://www.lofar.org/ 

^ http:/ /www. haystack. mit.edu/ast/arrays/mwa/ 
^ http:/ /astro.berkeley.edu/^dbacker/eor/ 



of sky provides a direct probe of the evolution in the neu- 
tral fraction of the IGM, and is referred to as the global 
step (|Shaver et al.lll999l : iGnedin fc Shaveij|2004l : iFurlanettd 
20061). Observation of individual Hll regions in the image 
regime will prob e quasar physics as wel l as the evolution of 
the neutral gas (IWvithe fc Loebl 12004 iKohler et all l2005l : 
iGefl fc Wvithell2008h . The most promising probe is the ob- 
servation of the power spectrum of fluctuations together 
with its evolution in redshift (see, e.g., iFurlanct to ct al] 
I2OO6I I. This observation would trace the evolution of neu- 
tral gas with redshift as well as the topology of the reion- 
isation process via the spatial depende nce of the statis- 
tics of redshifted 21-cm fluctuations (e.g ., Tozzi et al .ll200d: 
Furlanetto et al.ll2004l: iLoeb fc Zaldarria ga 2004: Iliev et al.l 



2006l : IWvithe fc Moralesll2007h . In this paper we focus on the 



statistical signature of reionisation in 21-cm power spectra. 

Due to theoretical uncertainties and a lack of observa- 
tional evidence the exact role quasars play in the reioni- 
sation process is uncertain. Although relatively short lived 
(tq ~ 10®~*yr) (see, e.g.. [Mar tini 2004, for a review), it is 
inferred t hat luminous quasars r adiate close to their Edding- 
ton limit (|Kollmeier et al.ll2006l ) such that the ionising pho- 
ton emission rate of luminous quasars can be~ 10"s"\The 
observed quasar population did not supply suffi cent UV pho- 
tons to reionise the Universe at z > 6 (e.g.. iFan et al.ii200ll : 
iDiikstra et al] |2004 iMeiksinI |2005| ) even though quasars 
dominate the production of ionising photons a.t z <. 3. How- 
ever, luminous quasars are known to exist at the edge of the 
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reionisation epoch (|Fan et al.ll2006l ). Moreover, these high- 
redsh ift quasars have a large clustering bias l|Shen et al.l 
I2OO7I ). implying that the effect of their ionising contribution 
on the statistical signature of the epoch of reionisation may 
be significant during the latter part of the pre-overlap era. 
The aim of this paper is to investigate the modification of 
the 21-cm power spectrum by quasars during the epoch of 
reionisation. Our results suggest that detailed simulations 
and reionisation models will need to consider quasar con- 
tribution to 21-cm power spectra when interpreting future 
observations. 

This paper is organised as follows: we begin in Sec- 
tion [2] by describing the inclusion of stellar reionisation 
within a semi-numerical scheme to produce realistic fiuctu- 
ations in the ionised gas distribution, here we also describe 
the method by which we populate our simulation field with 
quasars; in Section|3]we dicuss the choice of parameter space 
in which we explore the effect of quasars on the ionisation 
state of the IGM; Section |4] describes the statistical signa- 
tures used to characterise reionisation; our results, including 
their epoch dependence, are presented in Section [S] we dis- 
cuss our conclusions and their consequences for redshifted 
21-cm power spectrum measurements in Section [6] 

Throughout this paper we adopt the set of cosmolog- 
ical parameters determined by the Wilkinson M icrowave 
Anisotropy Probe ( WMAP) l|Komatsu et aLlbOOSl ) for a flat 
ACDM universe: fim = 0.27 (dark matter and baryons); 
Q.A = 0.73 (cosmological constant); S7b = 0.046 (baryons); 
h — 0.7 (Hubble constant); Us — 1 (primordial spectrum 
index); as = 0.8 (primordial spectrum normalisation). All 
distances are in comoving units unless stated otherwise. 



2 SEMI-NUMERICAL IONISATION MODEL 

There are two contributors to the reionisation of the IGM 
that we consider here: stars and quasars. In this section, 
we begin by summarising the features of our semi-numerical 
model for the reionisation of a three-dimensional volume of 
the IGM by galaxies fSection l2.1|l . Only a brief d escription is 
given here-we direct the reader to lGeil fc Wvit hc (2008) for 
further details. We then discuss our method of populating 
a sample volume with quasars using a Monte Carlo algo- 
rithm (Section 12. 2|) before describing the three-dimensional 
realisation of the ionisation state of the IGM including the 
ionising influence of both stars and quasars (Section 12. 3|l . 



2.1 Stellar ionisation 

We begin by simulating the linear matter overdensity field 
5{x, z) = Pni{x, z)/pm — 1 luslde a periodic, comoving, cubic 
region of volume = L^, by calculating the density contrast 
in Fourier space S(k,z) corresponding to a ACDM power 
spectrum (jEisenstoin & Hu 1999) linearly extrapolated to a 
specified redshift z. 

The semi-analytic model used to compute the relation 
between the local dark matter overdensity and the ionisa- 
tion state of the IGM is based on the mo del described by 
IWvithe fc L"oebl ^007) and ,Wvithe fc Morales (2007.). Any 
model for the reionisation of the IGM must describe the 
relation between the emission of ionising photons by stars 
in galaxies and the ionisation state of the intergalactic gas. 



This relation is non-trivial as it depends on various internal 
parameters which may vary with galaxy mass. These pa- 
rameters include the fraction of gas within galaxies that is 
converted into stars and accreting black holes, the spectrum 
of the ionising radiation and the escape fraction of ionising 
photons from the surrounding interstellar medium as well as 
the galactic halo and its immediate infall region (see iLoebl 
2006, for a review). This relation also depends on intergalac- 
tic physics. In overdense regions of the IGM, galaxies will be 
over-abundant because small-scale fluctuations need to be of 
lower amplitude to form a galaxy when embedded in a larger 
scale overdensity ||Mo fc White! 1 199^ . On the other hand, 
the increase in the recombination rate in overdense regions 
counteracts this galaxy bias. The process of reionisation also 
contains several layers of feedback: radiative feedback heats 
the IGM and results in the suppression of low-mass galaxy 
formation which delays the completion of reionisation by 
lowering the local star formation rate. However, this effect 
is counteracted in overdense regions by the biased formation 
of massive galaxies. 

The evolution of the ionisation fraction by mass Qsfi,R 
of a particular region of scale R with overdensity Sr (at 
redshift z) may be written as 



dQsf,,R 
dt 



0.76 



on / 

^ di 

rf-Fcol(5fl,^,2,Mmin)' 



dt 



aeCnH 



1+Sf. 



Diiz) 



^ Dl(Zobs 



(1 + zf Qs^,R, {1) 



where A'^ion is the number of photons entering the IGM per 
baryon in galaxies, as is the case-B recombination coeffi- 
cient and C is the clumping factor (which we assume, for 
simplicity, to be a constant value of 2). The production 
rate of ionising photons in neutral regions is assumed to 
be proportional to the collapsed fraction _Fcoi of mass in 
halos above the minimum threshold mass for star formation 
(Afmin), whereas the minimum halo mass in ionised regions is 
limited by the Jeans mass in an ionised IGM (Mion). We as- 
sume Mmin corresponds with a virial temperature of 10* K, 
representing the hydrogen cooling threshold, and Afion to 
correspond with a virial temperature of 10^ K, representing 
the mas s below which infal l from an ionised IGM is sup- 
pressed ijPiikstra et al.l[2004h . In a region of comoving radius 
R and mean overdensity 5{z) = SDi(z) / Di{zohs) (specified 
at redshift z instead of t he usual z = 0), we use the extended 
Press-Schechter model ijBond et al.l Il99lh to calculate the 
relevant collapsed fraction. 



Fcoi{Sr, R, z) = erfc 



5c - 3r{z) 
2(o^ 



gal 



(2) 



where erfc(a:;) is the complimentary error function, rj\ is the 
variance of the overdensity field smoothed on a scale R and 
cr|jji is the variance of the overdensity field smoothed on a 
scale -Rgai, corresponding to a meiss scale of Mmin or Mion 
(both evaluated at redshift z rather than at z = Q). 

Equation ([l]) may be integrated in time as a function 
of Sr. At a specified redshift this yields the filling fraction 
of ionised regions within the IGM on various scales J? as a 
function of overdensity. This model predicts the sum of as- 
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trophysical effects to be dominated by gaiaxy bias, and ttiat 
as a resuit, overdense regions are reionised first. Tliis ieads 
to tiie growth of H II regions via a pliase of percolation dur- 
ing wliicli individual Hll regions overlap around clustered 
sources in overdense regions of the Universe. We may also 
calculate the corresponding average 21-cm brightness tem- 
perature contrast of gas relative to th e cosmic microwave 
background (CMB) (|Madau et al.lll997t ). 

Sn{5n,R) = 26(l-Q.„,«)(l + 5b)(l-|j) 

(3) 



0.15 1 + z 



1/2 



mK. 



Here, the spin temperature Ts, CMB brightness temper- 
ature Tj, baryonic overdensity Sh (which we assume to 
be equal to the dark matter overdensity) and ionisation 
fraction Qsji,r are calculated for a comoving location x 
at redshift z. We assume ^ during the epoch o f 
reionisation (|Ciardi fc Madau|[2003l : iFurlanetto et al.l 120061 ) 
and ignore the enhancement of brightness temperature fluc- 
tuations due to peculiar velocities in overden se regions 
iBharadwai fc Alii l2005l : iBarkana fc Loebl |2005| ). Peculiar 
velocitie s were included in the semi-numerical model of 
[Mesing er fc Furlanettd (|2007l ). who found their effect to be 
small on scales ^ 10 Mpc. 

Throughout this paper we consider a model in which 
the mean IGM is fuUy reionised by z = 6 (|Fan et al.ll2006l : 
iGnedin fc FaiilbOOel : fwhite et all 120031 ). In this model we 
assume that star formation proceeds in halos above the hy- 
drogen cooling threshold in neutral regions of the IGM. In 
ionised regions of the IGM star formation is assumed to 
be suppressed by radiative feedback. Having determined the 
value of the ionised fraction Qsp.,R as a function of over- 
density &ii and smoothing scale R, we may now construct 
the ionisat ion field. We employ a s i milar filtering algorithm 
to that of iMesinger fc Furlanettd (|2007h to determine the 
mass-averaged ionisation state of each voxel within the sim- 
ulation cube. This is done by repeatedly filtering the linear 
overdensity field using a real-space spherical 'top hat' of ra- 
dius Rf, on scales in the range from L to L/N at logarithmic 
intervals of width ARt/Rt = 0.1. This procedure is per- 
formed in Fourier space for computational efficiency. For all 
filter scales the ionisation state of each voxel is determined 
using Qsji,R and deemed to be fully ionised if Qsfi,R ^ 1. AH 
voxels within a sphere of comoving radius R centred on these 
positions are fiagged and assigned Qsj{,r = 1, while the re- 
maining non-ionised voxels are assigned an ionised fraction 
of Qsji,R{ min' where -Rf,min = L/N corresponds to the small- 
est smoothing scale. A voxel forms part of an H II region if 
Qsr,r > 1 on any scale R. 



2.2 Quasar ionisation 

The ionising emissivity of the known population of quasars 
diminishes rapidly beyond z ^ 3. Furthermore, bright 
quasars are unlikely to c ontribute signicantly to the ionis- 
ing b ackgrou nd at z > 5 ( Shapir o et al.|[T994l : iHaiman et al.l 
|2001j : Wvithe fc Loebi r2003). More recent resuhs have con- 
strained the relative quasar contribution (at z ~ 6) to 
the ionising background t o be <^ 14 per cent (see, e.g., 
ISrbinovskv fc Wvithel2o"07l ). However, while the reionisation 



of cosmic hydrogen is believed to have been dominated by 
stars, the exact contribution made by quasars is poorly un- 
derstood . Moreover, the la rge clustering bias of high-redshift 
quasars (|Shen et al.ll2007l ) implies that even a small contri- 
bution to reionisation by quasars could have a significant 
effect on the statistics of 21-cm fluctuations. In this section, 
we describe the method by which we populate our simu- 
lation field with quasars and include their contribution to 
reionisation in our model. 



2.2.1 Quasar population 

In this paper, we model the quasar population as follows: 
all quasars are assumed to reside in dark matter halos and 
a given halo can host no more than one (active or dormant) 
quasar at a time; we assume the luminosity of a quasar dur- 
ing its active phase to be related to the mass of its host 
halo, with a resulting ionisation field that is isotropic. A 
mass-luminosity relationship is p lausible on both theoretical 
and empirical ground s (see, e.g.. iMartini fc WeinberdboOll : 
iKollmeier et al.ll200d ). however, the model we advocate here 
is phenomenologically motivated and not intended to encom- 
pass the physics of super-massive black hole growth. We do 
not consider the effects of anisotropic ionising sources. 



2.2.2 Space density of quasars 

If active quasars reside in a fraction /on of halos, then the 
comoving number density of active quasars 4>(m, z) in host 
halos with masses between (m, m + dm) at redshift z can 
be written as (j]{m, z)dm = /onWh(m, 2:)ijm, where nh(m ,z) 
is the Sheth-Tormen mass function (jSheth et al.l 1200 ll ) at 
redshift z. However, the effect of quasars on the history of 
reionisation is sensitive to the cumulative effect of the ion- 
ising contribution from quasars on the IGM, and so we are 
more concerned with the probability that a quasar has irra- 
diated its local volume at some time in its past. We therefore 
write 



(m, z)dm = nh(m, z)dm 



(4) 



as the comoving number density of quasars (active or dor- 
mant) in host halos with masses between (m, m -I- dm) at 
redshift z. We include fossil quasar-generated H II regions 
by assuming they remain ionised after the quasar has be- 
come dormant. The validity of this assumption is discussed 
in Section [2]3TT] 



2.2.3 Populating the simulation field 

The strong clustering of galaxies in overdense regions im- 
plies that these galaxies, and therefore quasars, should trace 
the higher density regions of the IGM. To study the impact 
high-redshift quasars have on the ionisation state of the IGM 
we populate our simulated field with quasars using a sim- 
ple case of a random walk Metropolis-Hastings Monte Carlo 
Markov chain (MCMC) algorithm, which we now describe 
in its applied context. 

The likelihood of a quasar residing at a random location 
is proportional to the number density of dark matter halos. 
For small values of large-scale overdensity 5, this number 
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Figure 1. Example quasar-populated field showing the target dis- 
tribution 7r(X) (grey scale) and projected quasar positions [closed 
circles) in 100^ Mpc^ X 0.4 Mpc simulation slice. HII regions sur- 
rounding quasars are shown as open circles of comoving radius i?q 
corresponding to their extent in a fully neutral IGM. 



density is proportional to 1 -)- b{M, z)S, where h{M, z) is the 
galaxy bias. 

We denote the Markov chain (MC) by X = 
{Xi, X2 , X„^ } , where X„ is the nth quasar placed in the 
simulation field, of mass M and at position x„. Given the 
current position of the MC is Xt, a move to Y is proposed 
by randomly perturbing the current position using a spheri- 
cally symmetric proposal distribution, q{x, y) = q{\x — y\). 
We accept the proposed move and set Xt — Y based on the 
acceptance probability 

a{Xt,Y)=mm{l,rM), (5) 
where is the Metropolis ratio, found by calculating 
n{Y)q{Y,Xt) _ n{Y) 



■K(Xt)q{Xt,Y) ■K{Xt)' 



(6) 



The second equality in Equation ((6)1 follows from the sym- 
metry of the proposal distribution q{x,y). The target dis- 
tribution n{X) is given by 



MX) 



l + b{M,z)SR{X) 



(7) 



1 + b{M, z)Sr 

, max 

where Sr denotes the overdensity field smoothed on a scal^B 
R = 5 Mpc. When b{M, z) = 1, 7r(X) is simply the smoothed 
dark matter overdensity at a sample position x normalised 
with respect to the maximum smoothed dark matter over- 



* As features in the overdensity field of scale ^ R disappear with 
larger smoothing scales, so do features in the target distribution 



niX). 



density in the realised volumjf|. We calculate the value for 
bias using the fitting formula of ISheth et al.l (|200ll ). 

Populating the density field in the manner described 
above for a single mass scale generates a sample set (the 
MC) with a probability density which (in the limit Nq ~ N^) 
converges on the desired target distribution more efficiently 
than using a transition probability dependent upon n{Y) 
alone. Since the overdensity field (and therefore the target 
distribution) is multi-modal, an MCMC with a peaked pro- 
posal distribution (relative to the spatial scale of the modes 
in the target distribution) can become bound about a lo- 
cal mode and fail to fully explore other modes of significant 
probability. We avoid this problem by using a uniform pro- 
posal distribution, ensuring that the random walk can reach 
all potential halo locations. This also removes the necessity 
of a burn-in period. We execute this algorithm for a range 
of quasar host halo masses above a minimum Mh,min. The 
minimum mass is chosen to be larger than the cosm ological 
Jeans mass in an ionised IGM (Mjjf| (|Gnedinll2000h . 

The large volume simulations required to investigate 
small frequency Fourier modes (i.e. large spatial scales) in 
the 21-cm signal at times approaching overlap do not in 
general resolve this Jeans scale. The simulation has a mass 
and volume resolution of Mcdi and Vcou respectively, where 
A/ccii = PrnKcii. A host halo of mass Mh < Mceii may or 
may not contain a quasar which is capable of fully ionising 
the simulation cell. This sub-cell ionisation issue is treated 
by determining the minimum halo mass required to ionise 
a simulation cell. Ma- If Mh > Ma the cell is fully ionised 
by the quasar, if Mh < Ma the cell is only partially ionised 
and the quasar contribution for that cell (added to the stellar 
contribution) is Qq — V/VccU- The model used to calculate 
Ma is outlined in Section 12.31 For example, using a 128"^ 
grid with a sidelength L = 300 Mpc does not resolve Mj 
at z = 7, however, a 5.1 x 10^" Mq halo housing a quasar 
would fully ionise a volume V ^ Vcoii. The maximum mass 
scale considered here is Mh.max ~ W^^ Mq, although the 
effective maximum mass scale is ~ 2.2 x 10^^ Mq (which 
is the maximum mass of a halo expected to be found in a 
volume of 300^ Mpc^' at z « 6). 



2.3 Inclusion of quasars in the semi-numerical 
scheme 

In a region of comoving radius R containing a quasar, the 
cumulative number of ionisations per baryon will be larger 
than predicted by Equation We incorporate the ionis- 
ing effect of a quasar by first computing the fraction of the 
IGM within a region of comoving radius R centred on the 
comoving position x that has already been reionised by stars 
[using Equation ((2)1], we then add an additional ionisation 



^ Note that n{X) can be negative for b{M, z) > 1. In such cases 
we set 7r(Jf) = 0. Furthermore, the normalisation constant for 
Tr{X) cancels in Equation ^ and therefore is not strictly required. 
We include it here to give tt{X) a probabilistic interpretation (i.e. 
0^7r{X)^l). 

^ We maintain the same mass bins for all simulations in order to 
reduce numerical noise in comparisons. 
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fraction equal to the quasar's contribution, 



Qq(x,M) = 



i?q(M) 



(8) 



where i?q(Af) is the comoving radius of the Hll region cen- 
tred on Xq that would have been generated by the quasar 
alone in a fully neutral and uniform IGM. In contrast to 
stellar ionisation the quasar contribution comes from a point 
source, and so the contribution to Q originates from a single 
voxel only (rather than all voxels within \x — Xq\). Follow- 
ing this addition we filter the ionisation field as described in 
Section O 

Consider a quasar located within a dark matter halo 
of mass Mh, radiating with a line-of-sight ionising photon 
emission rate of N.y photons per second. Assuming a uni- 
form IGM and neglecting recombinations, the proper radial 
extent of the quas ars Hll region obs e rved alon g the line of 
sight at time t is (|White et al.ll2003l : iHaiman fc Gen.. 200^ : 
iFan et alJliooel ) 



R = 3.2a;, 



-1/3 



1057 10'^ yr) 



1/3 



1+Z 

7.5 



pMpc, 



(9) 



where xni is the neutral fraction and tagc = t — -R(tagc)/c 
is the age of the quasar corresponding to the time when 
photons that reach R at time t were emitted. The value of 
R is subject to large uncertainties, including the quasar life- 
time and luminosity, as well as the dumpiness and ionisation 
state of the surrounding IGM. Because of the method used 
to include the ionising eff'ect of quasars at locations beyond 
the sphere that would be generated by the quasar alone in 
a fully neutral and uniform IGM, we put xni ~ 1 when 
calculating R. 

We assume that the mass of the central black hole (BH) 
scales as a power law with the circular velo city Vc of the 
host halo, Mbh oc v2 (|Wvithe fc Loebll2003l ). The circular 
velocity of a halo of mass Mh at redshift z varies as Vc{z) oc 
mI^^(1 + zfl"^ l|Barkana fc Loebl bOOll ). which leads to a 
power law relation between BH mass and halo mass, Mbh oc 
\ -\- zf'l'^ . If the BH powering this quasar shines at a 
fraction r\ of its Eddington luminosity then we have Lq oc 
riA'Ihh oc 'qM^^^{l + z)"^^^, w here 7 = 5 for simple models 
of feedback-limited accretion (jWvithe fc Loebl 120031) . Since 

Ln OC N-, 



jq II-,, using Equation ^ leads to (in proper units) 



R oc 



J7iagc(^) 



-,1/3 



M,:/^i+.) 



7/6-1 



(10) 



In models where BH gro wth is limited by feedb ack over the 
gas dynamical time (e.g.. IWvithe fc Loebl[20o3 l. the quasar 
lifetime t. 

2 S> 1. Equating t 
tion (|10|1 gives 



oq is proportional to H, and so tq 
with t, 



oc (1 + z)-^/2 for 
q and substituting into Equa- 



7?q(Mh,2) = Ron 



1/3 



Afh 



1012 M, 



5/9 



1 



7 



-2/3 



(11) 



Here we adopt a quasar l uminosity fraction of 77 = 0.1 
( Martini fc Weinberd l200lf) an d an index value 7 = 5 
i Wvithe fc Pa dmanab hanI 120061 ). Since the relative quasar 
contribution to reionisation is expected to be <^ 14 per cent 
(jSrbinovskv fc WvithdlioOTl ). we have chosen to investigate 
fiducial relative quasar ionisation contributions of 10 per 
cent. We therefore choose Ro and Mh min such that at over- 
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Figure 2. /on parameter solutions for simulations SI, S2 and S3. 



lap (zov = 6 in our model) the ionising contribution from a 
population of quasars alone relative to the ionising contri- 
bution of galaxies takes a value Qq/Q* = 0.1. 

We reiterate that the prescription outlined in this sec- 
tion is intended to provide a parameterised model, loosely 
based on underlying physical models, rather than a complete 
physical model of all the processes involved. The parameter- 
isation chosen provides a general set of models within which 
various different possible physical models are contained. The 
plausibility of this modelling can be investigated by deter- 
mining the values of various physical parameters for each of 
our simulation, which we denote by SI, S2 and S3 (in order 
of decreasing minimum halo mass). We do this for quasar 
duty cycle /on by noting that the fraction of active quasars 
with absolute bolometric magnitude A4b ^ —26 is 



/on — 



V'(Xe 



-26, 



ani,(m,z) 



(12) 



where M_26 is the host halo mass corresponding to a quasar 
with absolute bolometric magnitude A^b = —26. In order 
to solve for /on we note that for z ^ 1, 



tq ^ 9.78 X 10^ /l-'nn/''Von(l + zY^'I'^JT. 



(13) 



Equating R in Equation (|9} with i?q in Equation (|11|) en- 
ables Mh to be written as a function of /on and the model 
dependent J?o, which can then be substituted in Equa- 
tion (IT2I1. We use MM v. ^ -26) = 0.65 x 10"^ Mpc"^ at 
z = 6 llFan et al.ll2001^ and N-XM^i = -26) = 2 x lO^s"^ 
jWhite et al.ll200 3l). Figure [2l shows the solution for /on from 
Equation (|12p graphically for _Ro ~ 3, 4 and 5 proper Mpc. 
Table [1] gives the values of these parameters used in simula- 
tions S1-S3. We find duty cycles of ~ 10"^-4 x 10"^, which 
at z ~ 6 corresponds to tq ~ 10^^-4 x 10^ yr. These are com- 
parable to the Salpeter time fo r doubling of BH mass and 
to estimates of quasar lifetime l|Martinill2003 ). 
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Figure 3. Evolution of the upper-limit global ionised fraction 
Q due to the ionising contribution from the stellar component 
of galaxies only and quasars only. Quasar-only contributions are 
shown using 7 = 5/3 and a Sheth-Tormen mass function. The 
vertical grey line shows the redshift of overlap for this model, 



2.3.1 Recombination and fossil H II regions 

The prescription outlined above does not consider the ef- 
fect of recombinations within quasar-generated Hll regions. 
By investigating recombinatio ns in inhomogen e ous f ossil 
ionised regions around quasars, iFurlanetto et al.l (|2008l ) es- 
tablished that most of the Hll regions formed during hy- 
drogen reionisation remain highly ionised throughout the 
entire reionisation process. They found that the ionising 
background from galaxies inside the fossil Hll regions, to- 
gether with any residual low-level emission from the BH 
past its bright quasar phase, efficiently suppresses recom- 
binations during hydrogen reionisation. Thus, even though 
quasars are thought to be transient, the assumption of fully 
ionised relic bubbles yields a good approximation at the level 
of our model. 



3 THE QUASAR CONTRIBUTION TO 
REIONISATION 

We investigate the ionising contribution of quasars and the 
21-cm power spectrum as a function of redshift based on 
the parameterised model in Equation (|lip . together with 
the halo mass function and a minimum halo mass Mh.min- 
Table [1] gives the parameter values used in our simulations, 
which correspond to the plots in Figure [3] which shows the 
resulting contributions of stars and quasars to reionisation. 
Figure [3] shows that the assumption of a lower minimum 
halo mass allows quasars to contribute more significantly at 
earlier phases of reionisation. 

Figure |4] shows the relative differential ionised frac- 
tion due to quasars (Mh/Q*)dQq/dMh and the cumulative 
ionised fraction due to quasars Qq/Q*(< Afh) as a function 
of Mh for simulations S1-S3. The differential plots demon- 
strate that the bulk of the ionisation by quasars is achieved 
by quasars in lower-mass host halos, despite larger bubbles 
being generated by quasars in halos of higher mass. This is 
due to the greater population of lower-mass halos. As men- 



tioned above, the parameter sets used for these simulations 
have been chosen so as to produce a cumulative contribution 
to reionisation by quasars of 10 per cent. We also carried 
out simulations with a cumulative relative ionisation con- 
tribution by quasars of ~ 3 per cent. The results of these 
simulations (not presented) demonstrate that the effect on 
the resulting power spectra scale in proportion to quasar 
contribution. 



4 STATISTICAL SIGNATURES OF 
REIONISATION 

As mentioned in Section [TJ statistical observations of the 
epoch of reionisation promise to provide a wealth of infor- 
mation about the properties of neutral hydrogen at high 
redshift as well as some of the fundamental astrophysics 
behind the reionisation process and the first luminous ob- 
jects. While density perturbations in the matter distribution 
mediate fluctuations in the 21-cm signal both prior to and 
following reionisation, during the reionisation era the rela- 
tion between the 21-cm power spectrum and the underlying 
matter power spectrum is complex and, in its late stages, 
is dominated by the formation of large ionis ed "bubbles" 
l|Furlanetto et al.ll2004l : iMcQuinn et al.ll2006l ). These bub- 
bles of ionised hydrogen imprint features on the 21-cm power 
spectrum that reflect the luminosity and clustering of ion- 
ising sources responsible for reionisation. The 21-cm power 
spectrum can be used in a statistical test to distinguish can- 
didate reionisation models, as well as to constrain the history 
and morphology of reionisation (Barkana 200^). 

In this paper, we employ the dimensionless power spec- 
trum A2i{k,z) as the key statistical measure, which is the 
contribution to the variance of the redshifted 21-cm bright- 
ness temperature contrast 5Tb{z) per logarithmic interval 
in wave-number k = 2n/X. This measure is related to the 
dimensional form of the power spectrum P2i{k, z) by 



AUk,z) 



{2n 



1 -i 

ATvk'^P2i{k,z) 



(14) 



P2i(fc, z) is estimated by averaging over all m modes of the 
Fourier transform (T) of T{z) in a thin spherical shell in 
fc-space. 



^ m 

P2,ik,z)^{\f{k,z)f), = -j2mk, 



(15) 



Note that we must use caution when employing this sta- 
tistical measure since the spherical symmetry of the signal 
is broken by redshift evolution over certain ranges of red- 
shift. This issue has been discussed by a number of au- 
thors (s 



Morales fc HewittI |2004 iBarkana fc Loebl 



20051 : Im 

cQuinn et al.|[200^ ) and its effect on the sensitiv- 



ity of the 21-cm pow er spectrum has been investigated by 
iMcQuinn et al.l (|2006 ). Our numerical power spectrum mea- 
surements are subject to sample variance, arising from the 
finite number of independent modes counted in each fe-shell, 
which corresponds to the finite number of independent wave- 
lengths \ = 2-K /k that can fit into the simulated volum^^- 



^ In general, a measured power spectrum will also have a compo- 
nent of Poisson noise. Poisson noise is a discreteness effect that 
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Table 1. Parameter values used in our simulations. Simulation labels refer to the panel and numbered datam in Figure |3] Mx denotes 
a halo mass in units of 10-^ Mq. {N = 128, L = 300 Mpc) 



SimS 




A/h,min 


Ro (pMpc) 


/on 


iq (10^ yr) 


N^iz = 7) 


N^iz = 8) 


Afq(2 = 9) 


Afq(2= 10) 


SI 


10% 


6.6 Mil 


5.1 


0.04 


5.8 


14 











S2 


10% 


3.8A/11 


3.9 


0.005 


0.73 


108 


10 








S3 


10% 


1.6A/10 


3.0 


0.001 


0.15 


1294 


252 


34 


2 




Figure 4. Differential (left) and cumulative (right) ionisation fraction due to quasars Qq as a function of Afji for the simulation models 
considered. 



5 RESULTS 

We begin our analysis by presenting one example of the 
evolution of the IGM during the reionisation era with and 
without the ionising effect of quasars. The left-hand and 
central left-hand panels of Figure [5] show multi-phase ioni- 
sation maps for simulation S2 at z = 6.6, 7, 7.4 and 8, with 
and without the ionising influence of quasars respectively. 
The right-hand panel shows the corresponding power spec- 
tra and sample variance (faint lower curves). It is evident, 
both in the slices shown as well as in the corresponding aver- 
aged global ionisation fraction, that the simulated volumes 
with quasars are more ionised than those without at the 
same redshift. The resulting impact of the quasar-generated 
contribution to ionisation on the 21-cm power spectra is 
to decrease the overall power on all scales. The effect is 
epoch dependent since the accumulated quasar contribution 
to reionisation is redshift dependent (as shown in Figure |3]). 
The modification of the spectra arises due to the redistribu- 
tion of fluctuation power between wave-modes, originating 
from the spatial deformation of the stellar-generated bub- 
bles by the quasar emission. In addition to the modification 
of power, the addition of quasar ionisations also modifies the 
slope of the power spectra. This effect is most easily seen by 



is present due to imperfect sampling of the field, for example 
the finite number of galaxies in survey, or the gridding procedure 
used in creating inital conditions for Af-body simulations. Poisson 
noise is not present in our power spectra since we create our mass 
distributions using a fluid approximation. Rather than randomly 
place discrete massive particles in the simulation box in accor- 
dance with the required power spectrum of density fluctuations. 



calculating the relative amplitude between the power spec- 
tra without and with quasars, 

^21,QSOl"'' 

The resulting relative power for simulations SI, S2 and 
S3 are shown in the left-hand panels of Figures [6] and [T] 
Figure [6] shows the relative power at different redshifts as a 
function of wave- number. Figure [7] shows cuts of the data as 
a function of redshift at three different scales corresponding 
to wave-numbers k ~ 0.1, 0.4 and IMpc. The dependence 
on scale is discussed further in Section [S] 



5.1 The effect of quasars at constant neutral 
fraction 

In the examples above, the addition of quasar ionisation 
lowers the neutral fraction at a fixed redshift, and hence the 
comparison in Equation (|16p is not made at an equivalent 
stage of reionisation. Therefore, we next isolate the effect of 
quasars as a function of neutral fraction by removing the 
difference in overall power. This is achieved by boosting the 
stellar ionising efficiency in the stellar-only simulations so 
as to match the mass-weighted global neutral fraction of 
the simulation including quasars. The resulting "boosted" 
realisations for simulation S2 (to match xhi = 0.03, 0.16, 
0.28 and 0.43 in the examples) are shown in the central 
right-hand panels of FigureO and their corresponding power 
spectra and sample variance in the right-hand panels. 

The effect of quasars on the simulated 21-cm power 
spectra as a function of neutral fraction (and Mh,min) is 
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NO quasars INCLUDING quasars NO quasars BOOSTED 




0.1 1 

[Mpc-l] 



Figure 5. Multi-phase ionisation maps for simulation S2 (see Table[TJ a,t z = 6.6, 7, 7.4 and 8 {bottom to top) with (second column) and 
without {first column) the ionising influence of quasars, and with boosted ionisation fraction at xm = 0.03, 0.16, 0.28 and 0.43 (bottom 
to top). Each slice has a side length of 300 Mpc and is 2.3 Mpc deep. Fourth column: Spherically averaged three-dimensional 21-cm 
brightness temperature power spectra with (black) and without (red) quasars, and boosted (blue). Corresponding errors are shown below 
the spectra (faint). Fiducial sensitivity curves are shown for 1000 hr integrations using the MWA (dash) and the MWA5000 (dot-dash) 
with four 8 MHz processed sub-bands. 



again determined by calculating the relative power, 

a2 /, N ^il.noQSOeSiaJHl) , , 

A,el(fc,XHl) = — Ta 77 T' 

^21,QSo(«^'2;hi) 

where, as mentioned, n^qso has been calculated using 
the boosted simulations so as to match the mass-weighted 
neutral fraction of the simulations including quasars. 

The resulting relative power for simulations SI, S2 and 
S3 are shown in the right-hand panels of Figures |6] and [T] 



Figure |5] shows the relative power at different neutral frac- 
tions as a function of wave- number. Figure [7] shows cuts of 
the data as a function of neutral fraction at three differ- 
ent scales corresponding to wave-numbers k ~ 0.1, 0.4 and 
IMpc. The effect of the quasars is scale-dependent, indi- 
cating that the additional bias of the quasars makes fewer 
larger bubbles than a corresponding increase in stellar flux. 
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Figure 6. Comparisons between power spectra for simulations SI {top), S2 (middle) and S3 (bottom) with and without the inclusion 
of effects by quasars for 6.6 ^ 2: ^ 10 (left) and 0.03 ^ xm <^ 0.73 (right). The horizontal line marks unity which indicates equal power 
for both scenarios. The vertical lines corresponds to fc 0.1, 0.4, lMpc~^, or a spatial scale of 60, 15, 6Mpc respectively. 



6 DISCUSSION 

Our results indicate that a 10 per cent contribution to reion- 
isation from quasars leads to a damping of the 21-cm power 
spectrum by a factor of up to 2 at constant redshift, and by 
up to 30 per cent at constant neutral fraction shortly before 
the end of reionisation (Az ~ 1). The level of damping is 
both a function of scale and redshift. 

Firstly, in all simulations and at all scales, the effect 
on the relative power is more prominent at later epochs 
when there is a larger quasar population. At these later 
epochs, the cumulative quasar contribution to reionisation 
(from both newly formed and fossil H II regions) is greater 
(see Figure [3)1. Moreover, the rate at which the modification 
evolves is more rapid in models with more massive host ha- 
los, resulting from rapid evolution of the halo mass function. 

The scale dependence of the relative power can be at- 
tributed to the effects of halo bias. As discussed, our re- 
sults show that quasars increase the ionisation fraction at 



fixed redshift. This leads to a decrease in the amplitude of 
the power spectrum late in the reionisation era, since the 
growth of Hll regions removes signal at a rate that coun- 
teracts the increase in power resulting from the H II region- 
induced structure in the IGM. This is the "fall" of 21-cm 
fluctuations described by iLidz et al.l (2008). The effect is 
larger on small scales, since the addition of quasars increases 
the size of the Hll regions, owing to quasars being biased to 
overdense regions which reionise first. This removes more of 
the sources of small-scale power. The effect is less significant 
on large scales, where removal of the source of large-scale 
power due to the ionisation of gas is partly counteracted by 
its increase from larger Hll regions. 

In the case of comparison at constant neutral frac- 
tion the opposite is true. By increasing the contribution of 
star formation to match the neutral fraction, the fraction 
of small-scale power removed by ionisation is the same for 
stars and quasars. As a result A^^,! ~ 1 for large k. However, 
quasars result in fewer larger ionised regions as they have a 
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k= 1.0 
k=OA 
/c=0.1 



_1 I Ll. 



SI 



S2 



S3 



10 




0.19 



0.44 0.62 



0.75 



Figure 7. Comparisons between 21-cm power spectra with and without including the ionising of effect of quasars as a function of redshift 
{left) and neutral fraction (right) for simulations with Qq/Q^jzov = O-l- Three plots are shown for each simulation: fc 1 (dot-dashed), 
k ^ 0.4 (dashed) and k 0.1 Mpc"'^ (solid red). Instrumental sensitivity levels of the MWA (wide bar) and MWA5000 (narrow bar) for 
1000 hr integrations are shown for k O.lMpc"^ only. 



higher bias than the galaxies. Late in reionisation this ac- 
celerates the decline in power (during the "fall" ) and results 
in Ajgi > 1 for intermediate values of k. Thus, the level of 
modification of the shape of the 21-cm power spectrum can 
be related to the host halo masses of the quasars. 

Any radio interferometer is subject to instrumental 
noise and has limited sensitivity arising from the finite vol- 
ume of the observation. We consider observational parame- 
ters corresponding to the design specications of the MWA, 
and of a hypothetical follow-up to the MWA (which we 
term the MWA5000). We assume a continuous, circularly 
symmetric distribution of 512 antenna tiles for the MWA. 
This distribution has a constant antenna density core of ra- 
dius rc « 28 m and an inverse square radial antenna density 
profile for r > Vc (giving a maximum radius of 0.75 km). 
Each antenna tile contains 16 cross-dipoles to yield an ef- 
fective collecting area of = 16(AV4) (the area is capped 
for A > 2.1m). The MWA 5000 is assumed to follow the 
basic design of the MWA. The quantitative differences are 
that we assume the MWA5000 to have 5000 tiles within a 
diameter of 2 km, with a fiat antenna density core of ra- 
dius 40 m. In each case, we assume one field is observed for 
an integrated time of 1000 hr. Following the work of Mc- 
Quinn et al. (2006) we assume that foregrounds can be re- 
moved over 8 MHz sub-bands, within a bandpass of 32 MHz 
[foreground removal therefore imposes a minimum accessible 
wave-number of k^in ~ 0.04[(1 + z)/7.5]~^ Mpc"^]. 

We compute the sensitivity with which the effect of 



quasars on the shape of the 21-cm power s pectum could be 
detect ed fol lowing the procedure o utlined bv lMcQuinn et aT 
i2PQd) and iBowman et"ai] (|2006l ') (see also IWvithe et al 
2008). Written in terms of the cosmic wave- vector fc = 
fc|l + k±, the resulting error in the 21-cm power spectrum 
per mode is 



SP2i{k) 



D^AD 

Bto n^(JJ,u) 



+ P2l(fc), 



(18) 



where Tsys ~ 250[(l+2)/7]^'® K is the system temperature of 
the instrument, D[z) is the comoving distance to the point of 
emission at redshift z, AD is the comoving depth of the sur- 
vey volume corresponding to the bandwidth B,to is the total 
integration time, rib ((7, J^) is the number density of baselines 
that can observe the visibility U, where U = k±D/2TT and 
A is the observed wavelength. 

Although the observed 21-cm power spectrum is not 
spherically symmetric, it is symmetric about the line of 
sight. This makes it possible to calculate the overall power 
spectral sensitivity of the radio interferometer using the 
Fourier modes contained within an infinitesimal annulus 
around the line of sight of constant (k,6), where cos(0) = 
k ■ z/k (z is the unit vector pointing in the direction of the 
line of sight). The power spectral sensitivity over such an 
annulus is given by 



(jp{k,e) 



sP2i{k,e) 

^N^(k,9y 



(19) 
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where Nin{k, 9) denotes the number of observable modes in 
the annulus (only modes whose line-of-sight components fit 
within the observed bandpass are included). In terms of the 
fc-vector components k and 9, the number of independent 
Fourier modes within an annulus of radial width dk and 
angular width d9 is N^{k,e) = 2TTk^V sm{9)dkd9/{2nf, 
where V = AD{X^ /Ae) is the observed volume. Averag- 
ing (Jp{k, 9) over 9 gives the spherically averaged sensitivity 
to the 21-cm power spectrum ap{k), 



Ak,9)' 



(20) 



In order to find the sensitivity in terms of A21 we use Equa- 
tion (|14|) . which gives 



cr^2 (fc, z) = 



(2ny 



A-Kk'' 



E: 



-1/2 



(21) 



Figure [5] shows the fiducial spherically averaged sensitiv- 
ity curves for the MWA and MWA5000 (within bins of 
Ak — k/10) for integration times of 1000 hr using four 
8 MHz processed sub-bands. The corresponding sensitivities 
are also included in the relative power plots shown in Fig- 
ure[7]for k ~ 0.1 Mpc~^ only, where the instrumental sensi- 
tivities are near their maximum. 

Our results suggest that (for k ~ O.lMpc"^), using a 
constant redshift comparison, the effect of quasars on the 21- 
cm power spectrum is comparable to the precision achiev- 
able by the MWA and the MWA5000 at z ^ 7 for all three 
simulations, and at 7.5 <^ z <^ 8 in simulation S3. The sen- 
sitivity of the MWA5000 is comparable to the efi'ect up to a 
higher redshift of z ~ 9.4 in simulation S3. When making a 
comparison at constant neutral fraction, the effect is compa- 
rable to the observational precision of the MWA5000 only 
at xm <, 0.25 (power suppression) in all simulations, and 
0.4 < a^Hi < 0.65 (power amplification) in simulation S3. 
Thus, with respect to modelling the observable 21-cm power 
spectrum close to the end of reionisation, our results indi- 
cate that the potential contribution to ionisation by quasars 
will need to be considered. 



7 SUMMARY 

In this paper, we have assessed the effect of high-redshift 
quasars on the 21-cm power spectrum during the epoch of 
reionisation. Our approach has been to implement a semi- 
numerical scheme to calculate the three-dimensional struc- 
ture of ionised regions surrounding massive halos at high 
redshift. We have included the ionising influence of lumi- 
nous quasars by populating a simulated overdensity field 
with quasars using a Monte Carlo Markov Chain algorithm. 
Different parameterisations of the quasar luminosity-halo 
mass relationship have been used to analyse the relative 
effect on spherically averaged power in 21-cm emission be- 
tween simulations with and without the ionising influence 
of quasars. This comparison was carried out for simulations 
at the same redshift, ranging between 6.6 ^ z ^ 10, as 
well as at the same global neutral fraction, ranging between 
0.03 ^ a^Hi <; 0.75. Our results show that the cumulative 
ionising effect of quasars can suppress the 21-cm power spec- 
trum by a factor of up to 2 at constant redshift, and by up 



to 30 per cent at constant neutral fraction shortly before the 
end {Az ~ 1) of reionisation, both as a function of scale and 
redshift. 

The possible contribution of quasars to reionisation will 
complicate efforts to interpret observed 21-cm power spec- 
tra. For example, quasar ionisation results in modification 
of the amplitude of the power spectra at fixed neutral frac- 
tion. If quasars are ignored, this will result in biased infer- 
ences regard ing galaxy form ation (e.g., using a method as 
described in lBarkanal l2008l ). Similarly, the modification of 
the slope will lead to incorrect inferences regarding bubble 
size and host mass based on power spectral analyses that do 
not allow for quasars. Thus, modelling of the 21-cm power 
spectrum is likely to be degenerate with the quasar popu- 
lation, and so will require input from observations of the 
high-redshift luminosity functions of stars and/or quasars. 
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